Non-Dissipative Saturation of the Magnetorotational Instability in Thin Disks 
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A new non-dissipative mechanism is proposed for the saturation of the axisymmetric magnetoro- 
tational (MRI) instability in thin Keplerian disks that are subject to an axial magnetic field. That 
mechanism relies on the energy transfer from the MRI to stable magnetosonic (MS) waves. Such 
mode interaction is enabled due to the vertical stratification of the disk that results in the discretiza- 
tion of its MRI spectrum, as well as by applying the appropriate boundary conditions. A second 
order Dufiing-like amplitude equation for the initially unstable MRI modes is derived. The solutions 
of that equation exhibit bursty nonlinear oscillations with a constant amplitude that signifies the 
saturation level of the MRI. Those results are verified by a direct numerical solution of the full 
nonlinear reduced set of thin disk magnetohydrodynamics equations. 

PACS numbers: 47.65.Cb, 43.35.Fj, 62.60.+V 



Introduction - The destabilizing effect of an axial mag- 
netic field on Couette flow has been discovered half a cen- 
tury ago [H, 0. However, the importance of that phe- 
nomenon to astrophysics has been recognized only three 
decades later by Balbus and Hawley [l], [l]- That mecha- 
nism, termed the magnetorotational instability (MRI), is 
indeed considered by many researchers as the main can- 
didate to hold the key to solving the problem of angu- 
lar momentum transfer in accretion disks, and has been 
thoroughly investigated through linear analysis as well as 
nonlinear magnetohydrodynamic (MHD) simulations un- 
der a wide range of conditions and applications. Because 
of its importance in accretion disks physics the under- 
standing of the MRI's saturation mechanisms and level 
is of utmost significance. Thus, Knobloch & Julien [HI 
have described analytically the nonlinear saturation of 
the MRI in a straight infinite vertically uniform channel 
with solid boundaries. Such configuration is character- 
istic of laboratory experiments rather than astrophysical 
conditions. Knobloch & Julien considered a developed 
stage of the MRI, far from its threshold and showed that 
as the presence of the solid boundaries supports radial 
pressure gradients, the latter act together with the vis- 
cous as well as Ohmic dissipation in order to modify the 
rotation shear that feeds the instability and thus satu- 
rating it. In a complementary work, Umurhan ct al. Q 
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performed a weakly nonlinear analysis of the MRI close 
to marginality in configurations similar to Q and showed 
that the MRI saturates due to dissipative effects to levels 
that scale with the square root of the magnetic Prandtl 
number. In the current letter a novel mechanism for the 
saturation of the MRI is proposed, that differs from the 
above two works in the following important two aspects: 
1. A true thin disk is considered that is characterized 
by vertically localized stratified mass density, subject to 
radiation boundary conditions, and 2. The proposed 
mechanism is non-dissipative. Indeed, it is shown that 
the nonlinear forcing of magnetosonic (MS) waves by the 
MRI results in the saturation of the latter, and leads to 
bursty nonlinear oscillation of its amplitude. 
The Physical Model - The thin disk asymptotic expan- 
sion procedure is applied to the MHD equations in or- 
der to study the weakly nonlinear evolution of the MRI. 
The underlying physical property of the system is the 
supersonic nature of the Keplerian rotation whose Mach 
number is proportional to 1/e, where e is the ratio of the 
disk's semi-thickness to its characteristic radius. Conse- 
quently, both steady-state as well as the perturbed vari- 
ables are scaled with well defined powers of the small 
parameter e. Such procedure has been employed in nu- 
merous studies of thin disk dynamics, [7j-[ll|. and has in 
particular been proven efficient in the realistic analysis 
of the discrete MRI spectrum in true thin disk geometry, 
[l^ . Thus, the steady-state configuration is character- 
ized by a Keplerian rotation r2(r) = r^'^/^ where ri(r) 
and r are normalized by the value of the rotation fre- 
quency at some radius i?o, and Rq, respectively, as well 
as by an axial magnetic field Bz{r) that is an arbitrary 
function of r. The axial steady-state structure of the 
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disk is determined from the force balance between the 
thermal pressure and the gravitational pull of the central 
object. Thus, assuming axially isothermal configuration 
results in the following normalized number density pro- 
file: n(r,C) N{r)j:{ri), where £(77) = e''''/^^ N{r) is 
an arbitrary function of r, rj = C,/H{r), = z/e is the 
stretched axial coordinate, and H{r) is the semi thickness 
of the disk. The latter [or alternatively the temperature 
profile T(r)] is an arbitrary function of r. 

As the axial variations of the perturbations are as- 
sumed to occur on a much smaller length scales than the 
corresponding radial changes, to lowest order the MHD 
equations arc given by the following set of reduced non- 
linear equations that depend parametrically on the radial 
coordinate (see [12[ for detailed derivation): 
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where time is scaled with the local inverse rotation fre- 
quency, cr, {vr,ve,Vz) and {br,be,bz) are the perturbed 
number density (scaled by N{r)) and the components 
of the perturbed velocity (scaled by the sound velocity, 
Cs(r)), and magnetic field (scaled by the steady-state ax- 
ial magnetic field, Bz{r), respectively. In addition, the 
local plasma beta is given by /3(r) = P^N {r)Cl{r) / Bl{r) 
where /3o is the beta value at i?o, and Cs(r) is propor- 
tional to the square root of the temperature. 
Small Perturbations - Linearizing the system of equa- 
tions (HI)-® the latter decouples into two sub-systems: 
the first one is obtained from eqs. and describes 

the evolution of two Alfven-Coriolis waves one of which 
is responsible for the MRI. The second sub-system is ob- 
tained from eqs. (|5|)-(|6l) and describes the MS modes. 
Liverts & Mond [13[ obtained a WKB solution for the 
MRI eigenvalues and eigcnfunctions for Gaussian strati- 
fied disks. However, modifying the steady state number 



density axial profile to E(ry) = sech?r] enables the ana- 
lytical solution of both sub-systems. Assuming that the 
perturbations evolve as e''*'*, the eigcnfunctions of the 
Alfven-Coriolis sub-system can be expressed in terms of 
the Legendre polynomials which leads to the following 
dispersion relation [l2| : 



AX^p^ = 0, (7) 



where /3^^ = k{k + l)/3,fc = 1,2,.... As can be in- 
ferred from eq. (7), due to the axial stratification of the 
steady-state configuration the unstable modes are quan- 
tized with the axial number k which is equivalent to the 
axial wave number for the axially uniform case. Of par- 
ticular interest is the fact that the number of unstable 
MRI modes increases with f3, the threshold for exciting 
k unstable modes being /3^^. Figure 1 depicts the emer- 
gence of more and more unstable modes as (3 is increased. 
It is further noticed that each point on the /3-axis with 
f3 = Pcr^k ~ li2,..., serves as a bifurcation point for 
two modes, namely an unstable MRI one with +7 and a 
stable mode with —7. Consequently, at the bifurcation 
points the eigenvalue A is zero with multiplicity 2. This 
fact will turn out to be of great significance in the weakly 
nonlinear analysis to be unfolded in the next sections. 
Turning to the MS sub-system, its spectrum is stable 
and continuous. The eigenfunctions may be expressed as 
linear combinations of the following pair of linearly inde- 
pendent functions, Q: /± = [( i-o/ (i-hO]^^/^(V±o, 

where <^ = tanhrj, and fi ~ y/l — X^. It can be easily 
proven that solutions that vanish at ^ = ±1 exist only 
for > 0, which indeed renders the MS modes stable. 
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FIG. 1: The bifurcation diagram for MRI modes for k=l,2,3 
as a function of p. Each MRI branch that is characterized by 
a mode number k is accompanied by a stable branch that is 
symmetric about the vertical axis. 

Weakly Nonlinear Analysis - The focus is put now 
on the portion of the /3-axis that is only slightly above 
PIj. = 2/3. As can be seen from Fig. 1, which depicts 
the solution of eq. (7) for the first MRI modes, in that 
range of values there is only one such mode, namely, for 
fc = 1. Defining the control parameter (5 <C 1 in the 
following way: f3 = fS^ + S, the growth rate of the single 
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unstable MRI mode may be calculated to lowest order 
in 5 from eq. (7) to be 7^ = 27(5/14. The expression 
for the eigenfunction of, say, the radial component of 
the perturbed magnetic field of that initially small MRI 
perturbation looks, therefore as follows, [l^ : 



br{r,,j,t)^bo{r)[Po{0-^Pi{OHt), 



(8) 



where a{t) = aoe'^^^^*, Pk is the Legendre polynomial of 
order fc, fc = 0, 1, ^ = tanhr], bo is an arbitrary function 
of r, and oq is the amplitude of the initial perturbation. 

The amplitude a of the growing MRI remains exponen- 
tial as long as the perturbation is small enough. However, 
as the perturbation grows exponentially, nonlinear effects 
kick in. Thus, it is clear from eqs. (5) and (6) that the 
pressure exerted by the perturbed magnetic field excites 
the MS modes which in turn alters the stability properties 
of the original MRI mode by affecting its axial convection 
as well as by modifying the Alfven velocity. As a result of 
that interaction of the growing MRI with the excited MS 
mode the amplitude of the former is no longer exponen- 
tial but becomes a different function of time. The aim 
of the weakly nonlinear analysis is therefore to derive an 
appropriate ordinary differential equation that describes 
the evolution in time of the amplitude a{t) of the single 
MRI mode. In order to do that it is recalled that the 
transition to instability (when /3 reaches f3l^ from below) 
occurs when the linearized system of equations has a dou- 
ble zero eigenvalue. As a result the sought after equation 
is expected to be of second order as opposed to first order 
equations that characterize systems that bifurcate to in- 
stability through a simple zero eigenvalue Thus, the 
appropriate amplitude equation is of the following form: 
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Equation (9) has been derived directly from eqs. (l)-(6) 
by suitable asymptotic analysis . However, for sake of 
brevity of the current presentation, a is determined below 
from the parameters that characterize the steady-state. 
Equation (9) has also been used by Arter [l^ in order 
to describe sawtooth oscillation in Tokamaks. The role 
of double eigenvalues has also been recognized by Stcfani 
and Gerbeth [l7| in order to study polarity reversals in 
mean-field dynamo models. 

In order to calculate a it is first noticed that cq. (9) 
shares two important features with the full set of the 
reduced nonlinear equations (l)-(6): 1. For 6 < the 
origin of the phase space (i.e., the a~ da/dt plane) is the 
only fixed point of the dynamical system described by 
eq. (9). This is a reflection of the fact that for (5 < the 
only steady-state solution of the set (l)-(6) is the orig- 
inal basic Keplerian rotation. 2. For 6 > the single 
fixed point at the origin turns into a saddle point while 
two additional fixed points emerge which are centers and 
are located in symmetrically opposite locations with re- 
spect to the origin. The latter occur due to nonlinear 



effects and are given by duc/dt = 0,a^ = 7^ /a. This 
kind of behaviour, once again, reflects the fact that once 
the MRI grows exponentially, the nonlinear terms give 
rise to a new stable steady-state which may be calcu- 
lated through eqs. (l)-(6). Having that in mind a clear 
strategy emerges for calculating the value of a , namely, 
by the amplitude of the non linear steady state solution 
of eqs. (l)-(6). 

The Nonlinear Steady-State - The first steady-state 
solution of the set of equations (l)-(6) is obviously given 
by setting all the physical variables to zero. That solution 
describes the original basic Keplerian flow without any 
perturbations. However, as the beta value protrudes into 
the unstable region in parameter space, an additional 
perturbed steady-state solution emerges. Finding that 
nonlinear steady state starts by realizing that ^^^(C) = 
w"''(C) = w"*(0 = Oj where the superscript ns denotes 
the nonlinear steady state solution. Consequently, the 
following single ordinary differential equation is derived 
for the radial component of the steady state magnetic 
field: 



(fb^ 

de 



3/3 



d 

d^ 



(6- 



db" 



0, 



-/3(l-e2)_(fo-)2 

(10) 

where, as written above, f3 = + 5 with = 2/3. An 
asymptotic solution of eq. (10) in the limit of small pos- 
itive 5 is now obtained by expanding h^" in the following 
power series in (or alternatively in powers series in 
7): 



6^(0 = + {V5ftizbf\i) 



(11) 



To lowest order in \/5 the solution is given by the right 
hand side of eq. (8), exactly as the first linear eigenfunc- 
tion. The coefficient is determined now by the solv- 
ability condition of the resulting inhomogeneous equation 
for br {£,)■ Applying for that purpose Fredholm's alter- 
native theorem yields: ^.i = -^^5/2. Recalling now that 
the center fixed point of the dynamical system described 
by eq. (9) is given by = 7^ /a, and that Oc is identi- 
fied with V^/ii, the value of a is readily computed to be 
a = 27V5(5. 

Results - Equation (9) is known in the literature as 
the undamped Duffing equation and its solutions may 
be expressed analytically in terms of the Jacobi ellip- 
tic functions. Of particular importance is the fact that 
that kind of Duffing equation is derivable from an Hamil- 
tonian and hence its solutions are bounded. This re- 
fiects indeed the saturation of the MRI. For initial condi- 
tions near the saddle point at the origin of the phase 
plane the period of the bursty oscillations is asymp- 
totically given by O — >■ ln{8/h)/2'^ as ft, -> 0, where 
h = a(a^ — 7^a^ -I- a'^/2)t=to/'i'^^- Consequently, a small 
change in the initial conditions may dramatically change 
the period of the nonlinear oscillations without affecting 
in a significant way their amplitude. 
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FIG. 2: Comparison between 6r(?? = 0,t) as obtained from 
eqs. (l)-(6) (full line), and the solution of eq. (9) (dashed 
line), both for 5 = 0.0052. 
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Conclusions - A novel non dissipative mechanism for 
the saturation of the MRI has been proposed by which 
the latter drives non resonantly MS waves in a bursty 
oscillatory manner. A second order Duffing equation is 
derived for the amplitude of the MRI. The cubic term in 
that equation reflects the saturation mechanism in which 
the driven MS waves modify the plasma beta in a peri- 
odic way. Since the system is only slightly supercritical, 
due to that modification the average modified plasma 
beta oscillates above and below the instability threshold 
and thus instigates the observed bursty nonlinear oscilla- 
tions. Analytical solutions of the proposed Duffing equa- 
tion coincide with the solutions obtained from the nu- 
merical simulations of the reduced nonlinear MHD thin 
disk equations. Furthermore, the same dynamical behav- 
ior repeats itself near all threshold points /J^^, while nu- 
merical simulations indicate that the bursty oscillations 
persist also far away from IJ]. Therefore, the new 
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saturation mechanism presented in the current work im- 
poses severe limitations on the efSciency of the MRI to 
directly generate significant levels of turbulence in thin 
accretion disks. 
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FIG. 3: The amplitude of the perturbed number density in 
the equatorial plane, i.e., a{r] = Q,t), as obtained from eqs. 
(l)-(6) for 5 = 0.0052. 



Figure 2 depicts a comparison between hr{rj = 0,i) as 
obtained from the solution of the full nonlinear reduced 
set of MHD equations (l)-(6) (full line), and the solution 
of eq. (9) (dashed line), both for 5 = 0.0052. For that 
value of 5 the growth rate is 7 = 0.1. Thus, even though 
(5 <C 1 the characteristic time for saturation is quite real- 
istic at ^ 7 orbital periods. The two solutions practically 
overlap during the first few periods of the nonlinear oscil- 
lations. The phase difference between the two solutions is 
noticeable only later on due to the accumulating effect of 
a slight difference in the period. The latter, as has been 
discussed above, depends on and is very sensitive to the 
initial conditions. The amplitude of the nonlinear oscil- 
lations, which signifies the saturation level of the MRI, 
is the same for both calculations for the entire calcula- 
tion time. The growth immediately after each minimum 
value corresponds to exponential growth with the corre- 
sponding 7 value, i.e., 0.1. As a demonstration of the 
mode interaction mechanism that results in the satura- 
tion of the MRI, Fig. 3 depicts the perturbation in the 
number density at the equatorial plane, as obtained from 
the solution of eqs. (l)-(6) for 5 = 0.0052. The periodic 
energy transfer from the unstable MRI fc = 1 mode to 
the corresponding driven MS wave is indeed apparent. 
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